Demonstration: Reading Data

This notebook explains the data output from the GEODYN Starlette-SLR run.

First, load the data using the reading functions.

%load_ext autoreload
%autoreload 2

import pandas as pd
import numpy as np

##################################
# INPUT PARAMETERS:
##################################
sat_file = 'starlette'
arc1 = '030914_2wk'
grav_id ='goco05s'
local_path = '/data/analysis/starlette_analysis/'
SAT_ID = 7501001

##################################
# PATH TO DENSITY MODEL RUN of Choice:
##################################
msis86_model = 'msis86'
path_to_msis86 = '/data/runs_geodyn/st/results/'+ grav_id+'_'+msis86_model+'/'


import sys
sys.path.insert(0, '/data/analysis/util_funcs/py_starlette/')
from a_ReadStarlette import ReadStarlette

AdjustedParams, Trajectory, Density, Resids = ReadStarlette(arc1,
                                                           sat_file,
                                                           grav_id,
                                                           local_path,
                                                           path_to_msis86)

The base file name for this arc is: st030914_2wk.goco05s

File exists: iieout
        /data/runs_geodyn/st/results/goco05s_msis86/IIEOUT/st030914_2wk.goco05s
File exists: ascii_xyz
        /data/runs_geodyn/st/results/goco05s_msis86/XYZ_TRAJ/st030914_2wk.goco05s
File exists: densityfil
        /data/runs_geodyn/st/results/goco05s_msis86/DENSITY/st030914_2wk.goco05s

 Loading data...

Parameter adjustment data loaded
Trajectory data loaded
Density data loaded
Residual data loaded

Adjusted Parameters

• The Adjusted Parameters are read in from the iieout file (Unit 6).

• These include the State Vector at the start of the ARC, drag coefficient (CD), the solar radiation coefficient (CR), and the adjusted empirical accelerations (GA).

  • the CD is controlled/output by the DRAG card.

  • the GA is ouput by the ACCEL9 card (the 9 parameter general acceleration)

• AdjustedParams is a dictionary in which the first dimension stores each geodyn run iteration:

AdjustedParams.keys()

dict_keys([1, 2, 3, 4, 5])

And the 2nd dimension stores all of the adjusted parameters:

AdjustedParams[5].keys()

dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])

Each adjusted parameter then has the associated statistics:

APRIORI  VALUE
PREVIOUS VALUE
CURRENT  VALUE
TOTAL DELTA
CURRENT DELTA
APRIORI SIGMA
CURRENT SIGMA

AdjustedParams[5]['0XPOS']['CURRENT_VALUE']

-1228610.639773413

In the case of CD and GA, there are multiple values for each, and an additional dictionary exists for each value

print(AdjustedParams[5]['0GA'].keys(),'\n')

print(np.squeeze(list(AdjustedParams[5]['0CD'].keys())))

dict_keys(['0GA 9P 11', '0GA 9P 12', '0GA 9P 21', '0GA 9P 22'])

[Timestamp('2003-09-14 08:00:00') Timestamp('2003-09-14 16:00:00')
 Timestamp('2003-09-15 00:00:00') Timestamp('2003-09-15 08:00:00')
 Timestamp('2003-09-15 16:00:00') Timestamp('2003-09-16 00:00:00')
 Timestamp('2003-09-16 08:00:00') Timestamp('2003-09-16 16:00:00')
 Timestamp('2003-09-17 00:00:00') Timestamp('2003-09-17 08:00:00')
 Timestamp('2003-09-17 16:00:00') Timestamp('2003-09-18 00:00:00')
 Timestamp('2003-09-18 08:00:00') Timestamp('2003-09-18 16:00:00')
 Timestamp('2003-09-19 00:00:00') Timestamp('2003-09-19 08:00:00')
 Timestamp('2003-09-19 16:00:00') Timestamp('2003-09-20 00:00:00')
 Timestamp('2003-09-20 08:00:00') Timestamp('2003-09-20 16:00:00')
 Timestamp('2003-09-21 00:00:00') Timestamp('2003-09-21 08:00:00')
 Timestamp('2003-09-21 16:00:00') Timestamp('2003-09-22 00:00:00')
 Timestamp('2003-09-22 08:00:00') Timestamp('2003-09-22 16:00:00')
 Timestamp('2003-09-23 00:00:00') Timestamp('2003-09-23 08:00:00')
 Timestamp('2003-09-23 16:00:00') Timestamp('2003-09-24 00:00:00')
 Timestamp('2003-09-24 08:00:00') Timestamp('2003-09-24 16:00:00')
 Timestamp('2003-09-25 00:00:00') Timestamp('2003-09-25 08:00:00')
 Timestamp('2003-09-25 16:00:00') Timestamp('2003-09-26 00:00:00')
 Timestamp('2003-09-26 08:00:00') Timestamp('2003-09-26 16:00:00')
 Timestamp('2003-09-27 00:00:00') Timestamp('2003-09-27 08:00:00')
 Timestamp('2003-09-27 16:00:00') Timestamp('2003-09-28 00:00:00')]

AdjustedParams[5]['0GA']['0GA 9P 11']

{'APRIORI_VALUE': 0.0,
 'PREVIOUS_VALUE': 5.536419290656419e-11,
 'CURRENT_VALUE': 5.546253541325767e-11,
 'TOTAL_DELTA': -5.54625354133e-11,
 'CURRENT_DELTA': -9.83425066935e-14,
 'APRIORI_SIGMA': 1.0,
 'CURRENT_SIGMA': 1.6930522e-11}

AdjustedParams.keys()

dict_keys([1, 2, 3, 4, 5])

And the 2nd dimension stores all of the adjusted parameters:

AdjustedParams[5].keys()

dict_keys(['0CD', '0GA', '0GA_t2', '0XPOS', '0YPOS', '0ZPOS', '0XVEL', '0YVEL', '0ZVEL'])

import plotly.graph_objects as go
import numpy as np

Cd_plot = []

labels = list(AdjustedParams[5]['0CD'].keys())

for i in AdjustedParams.keys():

    Cd_plot.append(float(AdjustedParams[i]['0CD'][labels[0]]['CURRENT_VALUE']))

fig = go.Figure(data=[go.Scatter(x=list(AdjustedParams.keys()),
                                   y=Cd_plot,
                               mode='markers',
                                    marker=dict(
                                    size=10,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Drag Coefficient vs Iteration",
    yaxis_title=r"$C_d \,\,\text{(Unitless)}$",
    xaxis_title="Iteration Number",
    )

fig.show()

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline
from plotly.subplots import make_subplots


last_iter = 5
which_stat = 'CURRENT_VALUE' # 2 is the current val

labels = list(AdjustedParams[5]['0CD'].keys())


fig = make_subplots(
    rows=1, cols=1,
    subplot_titles=("Time Dependent Drag Coefficient",
                    ))
val_list = []
for i in AdjustedParams[last_iter]['0CD'].keys():
    val_list.append(AdjustedParams[last_iter]['0CD'][i][which_stat])
fig.add_trace(go.Scatter(x=labels,
                                 y=val_list,
                                  name= 'MSIS 86',
                               mode='markers',
                                    marker=dict(
                                    size=10,),
                                ), row=1, col=1,
                                   )


fig.update_yaxes( title=r"$C_D (Unitless)$",exponentformat= 'power',row=1, col=1)
fig.update_xaxes( title="Iteration Number", row=1, col=1)

# fig.update_yaxes( yaxis_title=r"$\frac{m}{s^2}$", xaxis_title="Iteration Number",exponentformat= 'power',row=2, col=1)

# fig.update_layout(
#     title="Along Track Accleration Amplitude vs Iteration",    )
# fig.update_layout(
#     yaxis_title=r"$\frac{m}{s^2}$",
#     xaxis_title="Iteration Number" )


iplot(fig)

Trajectory Data

• The trajectory data can come from multiple sources. The official output should come from the binary orbit file, but that is not working as of January 2020. Instead, we use the orbit file printout (ORBTVU card) to printout the trajectory data to an ascii file (Unit 10 for keplerian and Unit 08 for cartesian). These are then placed in the KEP_TRAJ and XYZ_TRAJ data directories before the temporary files are purged.

• The Trajectory structure is a dictionary in which the first dimension represents the different satellites, and the second dimension is a Pandas DataFrame containing the trajectory data from the ORBTVU card in which the first dimension stores each geodyn run iteration:

Indexing the first dimension will show the satellites in the GEODYN run that were saved. They are saved according to their SAT_ID. Since the Starlette-SLR runs only contain the Starlette satellite, this is the only one saved.

Trajectory.keys()

dict_keys([7501001])

The trajectory data can then be indexed by SAT_ID to show the pandas dataframe containing the data. The stored variables can be found with the .columns method.

print(Trajectory[7501001].columns)

Index(['YYMMDD', 'HHMM', 'SECONDS', 'X', 'Y', 'Z', 'XDOT', 'YDOT', 'ZDOT',
       'LAT', 'LONG', 'HEIGHT', 'timeHHMM', 'year', 'month', 'day', 'hours',
       'minutes', 'secs', 'Date'],
      dtype='object')

Trajectory[7501001].head()

	YYMMDD	HHMM	SECONDS	X	Y	Z	XDOT	YDOT	ZDOT	LAT	LONG	HEIGHT	timeHHMM	year	month	day	hours	minutes	secs	Date
2	30914	105	55.205719	-3902733.45	5158963.35	-3130061.00	-5675.94599	-1433.45008	4725.70967	-25.954786	118.024110	812276.80	0105	2003	09	14	01	05	55	2003-09-14 01:05:55
3	30914	106	15.005953	-4014288.94	5129496.49	-3035836.70	-5591.75463	-1542.86306	4791.45070	-25.120244	118.880387	811968.94	0106	2003	09	14	01	06	15	2003-09-14 01:06:15
4	30914	106	47.006428	-4190987.23	5077312.08	-2880862.22	-5450.70994	-1718.31837	4893.42231	-23.758977	120.237767	811607.96	0106	2003	09	14	01	06	47	2003-09-14 01:06:47
5	30914	107	13.806875	-4335433.07	5029307.46	-2748618.90	-5327.95616	-1863.82658	4974.68198	-22.607839	121.350810	811436.43	0107	2003	09	14	01	07	13	2003-09-14 01:07:13
6	30914	107	47.207484	-4510756.72	4964053.31	-2580843.47	-5169.22701	-2043.14535	5070.56125	-21.160247	122.709738	811390.75	0107	2003	09	14	01	07	47	2003-09-14 01:07:47

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline


#                     )
x = Trajectory[7501001]['Z'].astype(float)[:]
y = Trajectory[7501001]['X'].astype(float)[:]
z = Trajectory[7501001]['Y'].astype(float)[:]
fig = go.Figure(data=[go.Scatter3d(x=x,
                                   y=y,
                                   z=z,
                               mode='markers',
                                    marker=dict(
                                    size=2,
                                    opacity=0.7),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="3D Cartesian Orbit",
    xaxis_title="X Coordinate",
    yaxis_title="Y Coordinate",
#     zaxis_title="z Coordinate",
    )

# fig.show()
iplot(fig)

Density Data

• GEODYN does not actually have a built in method for printing out the density. As such, the density has been printd by modifying the DRAG.f90 subroutine in the IIE/MODS directory. The modification consists of printing the density and trajectory data at each timestep of the GEODYN run.

• The density data is stored as an ascii file named fort.99, which is then stored to the DENSITYdirectory before the temporary files are purged.

• The density data saves the longitutde, latitude, altitude coordinates, and cartesian coordinates.

Density.columns

Index(['Elapsed Secs', 'YYMMDD', 'HHMMSS', 'Lat', 'Lon', 'Height (meters)',
       'rho (kg/m**3)', 'delta_rho (kg/m**3/m)', 'X', 'Y', 'Z', 'XDOT', 'YDOT',
       'ZDOT', 'timeHHMMSS', 'year', 'month', 'day', 'hours', 'minutes',
       'secs', 'Date'],
      dtype='object')

Density.head()

	Elapsed Secs	YYMMDD	HHMMSS	Lat	Lon	Height (meters)	rho (kg/m**3)	delta_rho (kg/m**3/m)	X	Y	...	YDOT	ZDOT	timeHHMMSS	year	month	day	hours	minutes	secs	Date
0	210.0	30914	434	45.6973	278.5740	1083058.594	1.623031e-15	0.0	182992.974	-5215697.130	...	-1841.620957	-1912.851416	000434	2003.0	9.0	14.0	0.0	4.0	34.0	2003-09-14 00:04:34
1	225.0	30914	449	45.3666	279.6052	1084271.860	1.606673e-15	0.0	284238.712	-5242756.035	...	-1766.174236	-1989.534319	000449	2003.0	9.0	14.0	0.0	4.0	49.0	2003-09-14 00:04:49
2	240.0	30914	504	45.0256	280.6232	1085455.504	1.590705e-15	0.0	385422.960	-5268680.634	...	-1690.382976	-2065.747643	000504	2003.0	9.0	14.0	0.0	5.0	4.0	2003-09-14 00:05:04
3	255.0	30914	519	44.6746	281.6279	1086609.356	1.575121e-15	0.0	486523.867	-5293465.880	...	-1614.263040	-2141.476196	000519	2003.0	9.0	14.0	0.0	5.0	19.0	2003-09-14 00:05:19
4	270.0	30914	534	44.3139	282.6193	1087733.252	1.559916e-15	0.0	587519.621	-5317106.961	...	-1537.830323	-2216.704901	000534	2003.0	9.0	14.0	0.0	5.0	34.0	2003-09-14 00:05:34

5 rows × 22 columns

data_nums = 1000


import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density['Date'][:data_nums],
                                 y=Density['rho (kg/m**3)'][:data_nums],
                                    mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Density along Starlette Orbit (for 4 minutes of first day)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)

data_nums = 10
import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=Density['Date'][::data_nums],
                                 y=Density['rho (kg/m**3)'][::data_nums],
                                    mode='markers',
                                    marker=dict(
                                    size=2,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Density along Starlette Orbit (Every 10th datapoint)",
    yaxis_title=r'$\frac{kg}{m^3}$',
    xaxis_title="Date",
    )
fig.update_yaxes(type="log", exponentformat= 'power',)

iplot(fig)

Residual Data

• Similar to the trajectory data, the official output of the residuals should come from the binary residual file (Unit 19), but that is not working as of January 2020. Instead, we read in the residuals that are printed in ascii formats to the iieout file using the OBSVU card.

• The residuals are saved as a Pandas DataFrame.

Resids

	index	YYMMDD	HHMM	Sec-UTC-R	Observation	Residual	Ratio to sigma	Elev1	Elev2	OBS No.	...	timeHHMM	year	month	day	hours	minutes	secs	microsecs	Date	Ratio_to_sigma_fixed
0	0	30914	105	55.208831	933.052599576	-0.001484	-0.0015	58.485935	58.485839	1	...	0105	2003	09	14	01	05	55	208831	2003-09-14 01:05:55.208831	-0.0015
1	1	30914	106	15.009306	1005.251863666	-0.000437	-0.0004	51.260530	51.260437	2	...	0106	2003	09	14	01	06	15	009306	2003-09-14 01:06:15.009306	-0.0004
2	2	30914	106	47.010250	1145.739853454	-0.002657	-0.0027	41.550447	41.550361	3	...	0106	2003	09	14	01	06	47	010250	2003-09-14 01:06:47.010250	-0.0027
3	3	30914	107	13.811144	1279.798448021	-0.006932	-0.0069	35.007633	35.007553	4	...	0107	2003	09	14	01	07	13	811144	2003-09-14 01:07:13.811144	-0.0069
4	4	30914	107	47.212357	1461.073729804	-0.011274	-0.0113	28.408244	28.408171	5	...	0107	2003	09	14	01	07	47	212357	2003-09-14 01:07:47.212357	-0.0113
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3999	5688	30927	2306	16.259539	1260.651605954	0.006996	0.0070	53.959755	53.959671	4000	...	2306	2003	09	27	23	06	16	259539	2003-09-27 23:06:16.259539	0.0070
4000	5689	30927	2306	46.014220	1374.565632664	0.002611	0.0026	46.639968	46.639884	4001	...	2306	2003	09	27	23	06	46	014220	2003-09-27 23:06:46.014220	0.0026
4001	5690	30927	2307	14.118743	1496.430957455	-0.001994	-0.0020	40.669100	40.669018	4002	...	2307	2003	09	27	23	07	14	118743	2003-09-27 23:07:14.118743	-0.0020
4002	5691	30927	2307	46.624060	1650.422940823	-0.006786	-0.0068	34.755223	34.755145	4003	...	2307	2003	09	27	23	07	46	624060	2003-09-27 23:07:46.624060	-0.0068
4003	5692	30927	2308	5.727222	1745.985693706	-0.003580	-0.0036	31.700982	31.700906	4004	...	2308	2003	09	27	23	08	5.	27222	2003-09-27 23:08:05.027222	-0.0036

4004 rows × 21 columns

import plotly.graph_objects as go
import numpy as np
from plotly.offline import plot, iplot
%matplotlib inline

fig = go.Figure(data=[go.Scatter(x=pd.to_datetime(Resids['Date']),
                                 y=Resids['Residual'].values.astype(float),
                                    mode='markers',
                                    marker=dict(
                                    size=4,
                                    ),
                                   ),
                                   ],
                                   )

fig.update_layout(
    title="Observation Residuals from Final Iteration",
    yaxis_title="Residuals",
    xaxis_title="Date",
    )

iplot(fig)